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Abstract 

The modeling of tsunami generation is an essential phase in under- 
standing tsunamis. For tsunamis generated by underwater earthquakes, 
it involves the modeling of the sea bottom motion as well as the result- 
ing motion of the water above it. A comparison between various models 
for three-dimensional water motion, ranging from linear theory to fully 
nonlinear theory, is performed. It is found that for most events the linear 
theory is sufficient. However, in some cases, more sophisticated theories 
are needed. Moreover, it is shown that the passive approach in which the 
seafloor deformation is simply translated to the ocean surface is not always 
equivalent to the active approach in which the bottom motion is taken into 
account, even if the deformation is supposed to be instantaneous. 

1 Introduction 

Tsunami wave modeling is a challenging task. In particular, it is essential 
to understand the first minutes of a tsunami, its propagation and finally the 
resulting inundation and impact on structures. The focus of the present paper 
is on the generation process. There are different natural phenomena that can 
lead to a tsunami. For example, one can mention submarine mass failures, 
slides, volcanic eruptions, falls of asteroids, etc. We refer to the recent review 
on tsunami science |28| for a complete bibliography on the topic. The present 
work focuses on tsunami generation by earthquakes. 

Two steps in modeling are necessary for an accurate description of tsunami 
generation: a model for the earthquake fed by the various seismic parameters, 
and a model for the formation of surface gravity waves resulting from the defor- 
mation of the seafloor. In the absence of sophisticated source models, one often 
uses analytical solutions based on dislocation theory in an clastic half-space for 
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the seafloor displacement [55]. For tlie resulting water motion, the standard 
practice is to transfer the inferred seafloor displacement to the free surface of 
the ocean. In this paper, we will call this approach the passive generation ap- 
proach. This approach leads to a well-posed initial value problem with zero 
velocity. An open question for tsunami forecasting modelers is the validity of 
neglecting the initial velocity. In a recent note, Dutykh et al. [Sj used linear 
theory to show that indeed differences may exist between the standard passive 
generation and the active generation that takes into account the dynamics of 
seafloor displacement. The transient wave generation due to the coupling be- 
tween the seafloor motion and the free surface has been considered by a few 
authors only. One of the reasons is that it is commonly assumed that the source 
details are not important Ben-Menahem and Rosenman [T] calculated the 
two-dimensional radiation pattern from a moving source (linear theory). Tuck 
and Hwang [33] solved the linear long- wave equation in the presence of a moving 
bottom and a uniformly sloping beach. Hammack |15| generated waves exper- 
imentally by raising or lowering a box at one end of a channel. According to 
Synolakis and Bernard ^S], Houston and Garcia [17] were the first to use more 
geophysically realistic initial conditions. For obvious reasons, the quantitative 
differences in the distribution of seafloor displacement due to underwater earth- 
quakes compared with more conventional earthquakes are still poorly known. 
Villeneuve and Savage [35| derived model equations which combine the linear ef- 
fect of frequency dispersion and the nonlinear effect of amplitude dispersion, and 
included the effects of a moving bed. Todorovska and Trifunac [31| considered 
the generation of tsunamis by a slowly spreading uplift of the seafloor. 

In this paper, we mostly follow the standard passive generation approach. 
Several tsunami generation models and numerical methods suited for these mod- 
els are presented and compared. The focus of our work is on modelling the fluid 
motion. It is assumed that the seabed deformation satisfies all the necessary 
hypotheses required to apply Okada's solution. The main objective is to confirm 
or infirm the lack of importance of nonlinear effects and/or frequency dispersion 
in tsunami generation. This result may have implications in terms of compu- 
tational cost. The goal is to optimize the ratio between the complexity of the 
model and the accuracy of the results. Government agencies need to compute 
accurately tsunami propagation in real time in order to know where to evacuate 
people. Therefore any saving in computational time is crucial (see for example 
the code MOST used by the National Oceanic and Atmospheric Administra- 
tion in the US [2^ or the code TUNAMI developed by the Disaster Control 
Research Center in Japan). Liu and Liggett [24| already performed compar- 
isons between linear and nonlinear water waves but their study was restricted 

^In the pioneering paper |18l . Kajiura analyzed the apphcabihty of the passive approach 
using Green's functions. In the tsunami hterature, this approach is sometimes called the 
piston model of tsunami generation. 

■^As pointed out by Geist et al. [8], the 2004 Indian Ocean tsunami shed some doubts 
about this belief. The measurements from land based stations that use the Global Positioning 
System to track ground movements revealed that the fault continued to slip after it stopped 
releasing seismic energy. Even though this slip was relatively slow, it contributed to the 
tsunami and may explain the surprising tsunami heights. 
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to simple bottom deformations, namely the generation of transient waves by an 
upthrust of a rectangular block, and the nonlinear computations were restricted 
to two-dimensional flows. Bona et al. [3] assessed how well a model equation 
with weak nonlinearity and dispersion describes the propagation of surface wa- 
ter waves generated at one end of a long channel. In their experiments, they 
found that the inclusion of a dissipative term was more important than the in- 
clusion of nonlinearity, although the inclusion of nonlinearity was undoubtedly 
beneficial in describing the observations. The importance of dispersive effects 
in tsunami propagation is not directly addressed in the present paper. Indeed 
these effects cannot be measured without taking into account the duration (or 
distance) of tsunami propagation [32] . 

The paper is organized as follows. In Section 2, we review the equations 
that are commonly used for water-wave propagation, namely the fully nonlinear 
potential flow (FNPF) equations. Section 3 provides a description of the linear 
theory, with explicit expressions for the free-surface elevation and the velocities 
everywhere inside the fluid domain, both for active and passive generations. 
Section 4 is devoted to the nonlinear shallow water (NSW) equations and their 
numerical integration by a finite volume scheme. In Section 5 we briefly describe 
the boundary element numerical method used to integrate the FNPF equations. 
The following section (Section 6) is devoted to comparisons between the various 
models and a discussion on the results. The main conclusion is that linear 
theory is sufficient in general but that passive generation overestimates the 
initial transient waves in some cases. Finally directions for future research are 
outlined. 

2 Physical problem description 

In the whole paper, the vertical coordinate is denoted by z, while the two hori- 
zontal coordinates are denoted by x and y, respectively. The sea bottom defor- 
mation following an underwater earthquake is a complex phenomenon. This is 
why, for theoretical or experimental studies, researchers have often used simpli- 
fied bottom motions such as the vertical motion of a box. In order to determine 
the deformations of the sea bottom due to an earthquake, we use the analytical 
solution obtained for a dislocation in an elastic half-space [35]. This solution, 
which at present time is used by the majority of tsunami wave modelers to 
produce an initial condition for tsunami propagation simulations, provides an 
explicit expression of the bottom surface deformation that depends on a dozen 
of source parameters such as the dip angle 6, fault depth df, fault dimensions 
(length and width). Burger's vector D, Young's modulus, Poisson ratio, etc. 
Some of these parameters are shown in figure |T] More details can be found in 
[4] for example. A value of 90° for the dip angle corresponds to a vertical fault. 
Varying the fault slip |D| does not change the co-seismic deformation pattern, 
only its magnitude. The values of the parameters used in the present paper are 
given in Table [T] A typical dip-slip solution is shown in figure [2] (the angle (f> is 
equal to 0, while the rake angle 9 is equal to n/2). 
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Figure 1: Geometry of the source model (dip angle 5, depth o?/, length L, width 
W) and orientation of Burger's vector D (rake angle 9, angle (j) between the 
fault plane and Burger's vector). 



parameter value 

Dip angle 5 13^ 

Fault depth df, km 3 

Fault length L, km 6 

Fault width W, km 4 
Magnitude of Burger's vector |D|, m 1 

Young's modulus E, GPa 9.5 

Poisson ratio v 0.23 



Table 1: Typical parameter set for the source used to model the seafloor de- 
formation due to an earthquake in the present study. The dip angle, Young's 
modulus and Poisson ratio correspond roughly to those of the 2004 Sumatra 
event. The fault depth, length and width, as well as the magnitude of Burger's 
vector, have been reduced for computation purposes. 
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Figure 2: Typical seafloor deformation due to dip-slip faulting. The parameters 
are those of Table [T] The distances along the horizontal axes x and y are 
expressed in kilometers. 
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Figure 3: Initial net volume V (in km^) of the seafloor displacement as a function 
of the dip angle S (in °). All the other parameters, which are given in Table [TJ 
are kept constant. 



Let z = C,{x,y^t) denote the deformation of the sea bottom. Hammack 
and Segur |16| suggested that there are two main kinds of behaviour for the 
generated waves depending on whether the net volume V of the initial bottom 
surface deformation 

^ = y" Cix, y,0)dxdy 

is positive or notH A positive V is achieved for example for a "reverse fault", 
i.e. when the dip angle S satisfies < S < tt/2 or — tt < S < — 7r/2, as shown in 
figure [3l A negative V is achieved for a "normal fault", i.e. when the dip angle 
6 satisfies tt/2 < 5 < tt oi —tt/2 < S < 0. 

The conclusions of [T2| are based on the Korteweg-de Vrics (KdV) equation 
and were in part confirmed by their experiments. If V is positive, waves of stable 
form (solitons) evolve and are followed by a dispersive train of oscillatory waves, 
regardless of the exact structure of ({x, y, 0). li V is negative, and if the initial 
data is non-positive everywhere, no solitons evolve. But, if V is negative and 
there is a region of elevation in the initial data (which corresponds to a typical 
Okada solution for a normal fault), solitons can evolve and we have checked 

•^However it should be noted that the analysis of 1161 is restricted to one-dimensional uni- 
directional waves. We assume here that their conclusions can be extended to two-dimensional 
bi-directional waves. 
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Figure 4: Wave profiles at different times for the case of a normal fault {5 = 
167°). The seafloor deformation occurs instantaneously at t = 0. The water 
depth h{x,y) is assumed to be constant. 

this last result using the FNPF equations (see figure |4|). In this study we focus 
on the case where V is positive with a dip angle S equal to 13°, according to 
the seismic data of the 26 December 2004 Sumatra- Andaman event (see for 
example [12] )■ However, the sea bottom deformation often has an A^— shape, 
with subsidence on one side of the fault and uplift on the other side as shown 
in figure m In that case, one may expect the positive V behaviour on one side 
and the negative V behaviour on the other side. Recall that the experiments of 
Hammack and Segur [16] were performed in the presence of a vertical wall next 
to the moving bottom and their analysis was based on the uni-directional KdV 
wave equation. 

We now consider the fluid domain. A sketch is shown in figure O The fiuid 
domain fl is bounded above by the free surface and below by the rigid ocean 
floor. It is imboimdcd in the horizontal x— and y— directions. So, one can write 

n = R^x [-hix, v) + C(x, y, t), 77(.T, y, t)]. 

Before the earthquake the fluid is assumed to be at rest, thus the free surface 
and the solid boundary are defined by z = and z = —h{x, y), respectively. For 
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simplicity h{x, y) is assumed to be a constant. Of course, in real situations, this 
is never the case but for our purpose the bottom bathymetry is not important. 
Starting at time t = 0, the solid boundary moves in a prescribed manner which 
is given by 

z ^ -h + Cix,y,t), t>0. 

The deformation of the sea bottom is assumed to have all the necessary prop- 
erties needed to compute its Fourier transform in x, y and its Laplace transform 
in t. The resulting deformation of the free surface z — ri(x, y, t) is to be found 
as part of the solution. It is also assumed that the fluid is incompressible and 
the flow irrotational. The latter implies the existence of a velocity potential 
(f){x,y,z,t) which completely describes the flow. By definition of (j) the fluid 
velocity vector can be expressed as q = V0. Thus, the continuity equation 
becomes 

V-q = A(/. = 0, {x,y,z)en. (1) 

The potential (j){x,y,z,t) must satisfy the following kinematic boimdary condi- 
tions on the free surface and the solid boundary, respectively: 

= v{x,y,t), 

= -h + C{x,y,t). 

Further assuming the flow to be inviscid and neglecting surface tension ef- 
fects, one can write the dynamic condition to be satisfied on the free surface 
as 

^ + ^\V(f>\^ +gi] = 0, z = 7ix,y,t), (2) 
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where g is the acceleration due to gravity. The atmospheric pressure has been 
chosen as reference pressure. 

The equations are more transparent when written in dimcnsionless variables. 
However the choice of the reference lengths and speeds is subtle. Different 
choices lead to different models. Let the new independent variables be 

x = x/X, y = y/X, z — z/d, t — cot/X, 

where A is the horizontal scale of the motion and o? a typical water depth. The 
speed Co is the long wave speed based on the depth d (cq = V9^)- Let the new 
dependent variables be 

11=-, C=-, 9= — 70, 
a a agX 

where a is a characteristic wave amplitude. 

In dimcnsionless form, and after dropping the tildes, the equations become 

dl^ + ^ [d^ + WJ^^' (3) 



(6) 

where two dimcnsionless numbers have been introduced: 

£ = a/d, ^ = d/X. (7) 

For the propagation of tsunamis, both numbers e and /i are small. Indeed the 
satellite altimetry observations of the 2004 Boxing Day tsunami waves obtained 
by two satellites that passed over the Indian Ocean a couple of hours after the 
rupture process occurred gave an amplitude a of roughly 60 cm in the open 
ocean. The typical wavelength estimated from the width of the segments that 
experienced slip is between 160 and 240 km |22| . The water depth ranges from 4 
km towards the west of the rupture to 1 km towards the east. Therefore average 
values for e and fji in the open ocean are e « 2 x 10"'' and « 2 x 10^^. A 
more precise range for these two dimcnsionless numbers is 

1.5 X 10""* < £ < 6 X 10"^ 4 X 10^3 <^l< 2.5 x 10"^ (8) 

The water-wave problem, either in the form of an initial value problem (IVP) 
or in the form of a boundary value problem (BVP), is difhcult to solve because of 
the nonlincaritics in the boundary conditions and the unknown computational 
domain. 
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3 Linear theory 

First we perform the linearization of the above equations and boundary condi- 
tions. It is equivalent to taking the limit of ([21)-® as e ^ 0. The linearized 
problem can also be obtained by expanding the unknown functions as power 
series of the small parameter e. Collecting terms of the lowest order in e yields 
the linear approximation. For the sake of convenience, we now switch back to 
the physical variables. The linearized problem in dimensional variables reads 

A(j) = 0, {x,y,z) eM^ X [-h,0], (9) 

ijr = i^' ^ = 0' (10) 



dz dt 
dz dt ' 



z = ~h, (11) 

^+5?7 = 0, z = 0. (12) 

The bottom motion appears in equation (jlip . Combining equations pop and 
yields the single free-surface condition 

d^S d(h , , 

^+5^-0, z^O. (13) 

Most studies of tsunami generation assume that the initial free-surface de- 
formation is equal to the vertical displacement of the ocean bottom and take a 
zero velocity field as initial condition. The details of wave motion are completely 
neglected during the time that the source operates. While tsunami modelers of- 
ten justify this assumption by the fact that the earthquake rupture occurs very 
rapidly, there are some specific cases where the time scale and/or the horizontal 
extent of the bottom deformation may become an important factor. This was 
emphasized for example by Todorovska and Trifunac [31] and Todorovska et 
al. [30], who considered the generation of tsunamis by a slowly spreading uplift 
of the seafloor in order to explain some observations related to past tsunamis. 
However they did not use realistic source models. 

Our claim is that it is important to make a distinction between two mecha- 
nisms of generation: an active mechanism in which the bottom moves according 
to a given time law and a passive mechanism in which the seafloor deformation 
is simply translated to the free surface. Recently Dutykh et al. [5] showed 
that even in the case of an instantaneous seafloor deformation, there may be 
differences between these two generation processes. 

3.1 Active generation 

Since in this case the system is assumed to be at rest at < = 0, the initial 
condition simply is 

r7(x,2/,0) = 0. (14) 
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In fact, r]{x,y,t) = for all times t < and the same condition holds for the 
velocities. For t < 0, the water is at rest and the bottom motion is such that 
C{x,y,t) = for t < 0. 

The problem ((9|)- p3)) can be solved by using the method of integral trans- 
forms. We apply the Fourier transform in {x,y), 

m = fikj) = / f{x,y)e-'^''--+'yUxdy, 



with its inverse transform 

S-'[f] = fix,y) = I fik,e)e^^'^^+'yUkdi, 

1 

and the Laplace transform in time i, 

-t-oo 

)e-"* dt. 



m^s{s)= J 9{t)c 







For the combined Fourier and Laplace transforms, the following notation is 
introduced: 

+ 00 

^£[F(2;, y, t)] = F{k, £,s) = J e-'('=-+^2') dxdy J F(x, y, t)e-^* dt. 

After applying the transforms, equations ([5]), pT|) and become 

+ ^2)^ = 0, (15) 

'^'^{k,e,-h,s) ^ sCik,£,s), (16) 



dz 



s''(f>ik,£,0,s)+g^ikj,0,s)^0. (17) 
The transformed free-surface elevation can be obtained from (IT^: 



r]{k,i,s) = --(j){k,£,0,s). (18) 
9 

A general solution of equation (fT5|) is given by 

(j){k, £, z, s) = A{k, £, s) cosh(mz) -I- B{k, £, s) sinh(mz), (19) 



where m = \/k^ + W. The fimctions A{k, £, s) and B{k, £, s) can be easily found 
from the boundary conditions (fTB|) and ([T7]): 

A{k,£,s) = 9<{k,£,.s) 



B{k,£,s) 



cosh(m/i)[s2 -|- gim tanh(m/i)] ' 

s^~ak,£,s) 

m cosh(TOft,) [s^ + gm tanh(mft,)] 
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From now on, the notation 



Lo = \/ gm tanh(mft.) (20) 



will be used. Substituting the expressions for the functions A and B in ([T 
yields 



s 



2 



cj)(k, i, z, s) = —, — , ' ' TTT cosh(TOz) sinh(mz) . (21) 



The free-surface elevation becomes 



77(fc,£,s) 



cosli(m/i)(s2 -I- tj2-) 



Inverting the Laplace and Fourier transforms provides the general integral 
solution 

/^-t-zoo 



rj{x,y,t) = —- ——-— / V ' \' e^'dsdkdt (22) 
(27rj^ J J cosh[mn) 2m J + uu^ 



^ — lOC 



In some applications it is important to know not only the free-surface eleva- 
tion but also the velocity field inside the fluid domain. In the present study we 
consider seabed deformations with the structure 



C(x,y,t) :=Co(x,2;)r(i). (23) 

Mathematically we separate the time dependence from the spatial coordinates. 
There are two main reasons for doing this. First of all we want to be able to 
invert analytically the Laplace transform. The second reason is more fundamen- 
tal. In fact, dynamic source models are not easily available. Okada's solution, 
which was briefly described in the previous section, provides the static sea-bed 
deformation C,o{x,y). Hammack |15| considered two types of time histories: an 
exponential and a half-sine bed movements. Dutykh and Dias 0] considered two 
additional time histories: a linear and an instantaneous bed movements. We 
show below that taking an instantaneous seabed deformation (in that case the 
function T{t) is the Hcavisidc step function) is not equivalent to instantaneously 
transferring the seabed deformation to the ocean surfaccQ. 

In equation (|2ip. we obtained the Fourier-Laplace transform of the velocity 
potential 0(a;, z, t): 

— gsQ){k,£)T{s) f \ 

6(k, £, z, s) = '—, — , ' „ ^ cosh(TOz) sinh(mz) ) . (24) 

*In the framework of the Unearized shallow water equations, one can show that it is equiv- 
alent to take an instantaneous seabed deformation or to instantaneously transfer the seabed 
deformation to the ocean surface I33| . 



3.1 Active generation 



13 



Let us evaluate the velocity field at an arbitrary level z = Ph with — 1 < /? < 0. 
In the linear approximation the value /? = corresponds to the free surface 
while P = —I corresponds to the bottom. Below the horizontal velocities are 
denoted by u and the horizontal gradient {d/dx,d/dy) is denoted by V/j. The 
vertical velocity component is simply w. The Fourier transform parameters are 
denoted by k = (fc,£). 

Taking the Fourier and Laplace transforms of 

u{x,y,t;/3) = ^ h<l>{x,y, z,t)\^^p^ 

yields 

u{kJ,s;P) = -~i^{k,e,f3h,s)k 
. gsCoik,e)T{s) 



cosh{(3mh) smh{(3mh) ) k. 



cosh(m/i)(s2 + Lj2) \ gm 

Inverting the Fourier and Laplace transforms gives the general formula for the 
horizontal velocity vector: 

ig ff kCojkJ) cosh(m/3/i)e'(""+^^) 1 / sT{s)e^ 

]R2 fi—ioc 



^ /i + ZOO 

kCo(fc,^)sinh(m/3/i)e*('=^+^y) 1 f s^T{s)e^^ 



47r^ J J m cosh(m/i) 27ri 

]R2 fi — ioo 



ds dk. 



Next we determine the vertical component of the velocity w(x,y,t;f3). It 
is easy to obtain the Fourier-Laplace transform w{k,£, s; P) by differentiating 
(El: 



w{k,i,s]P) = — 



sgCo{k,e)T{s) 

z=l3h COsh(m/l)(s2 + tj2-) yg 



— cosh(/3TO/i) — m sinh(/3m/i) 



Inverting this transform yields 

w[x,y,t,p) ^^^jj ^^^^^f^) 2m J + lo^ 

]R2 fi—ioo 

^ /j,+'ioo 

j_ /■/■ msinh(/3mfe)Co(fc,£) 1 f sT{s)e/ ' 



'_^i(Kx+ty) I ds dk, 

47r^ J J cos\\{mh) 27ri J + ' 

for -1 < /3 < 0. 

In the case of an instantaneous seabed deformation, T{t) = H{t), where 
H{t) denotes the Heaviside step function. The resulting expressions for -q, u 
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and w (on the free surface), which are valid for i > 0, are 

- wfjj :X,k) <^^) 

u x,y,t;0) = / / , . , > c?k, (26) 

47r"^ J J cosh(m/i) w 

= -^// ^°'j,(,„,) ^^..^'^k. (27) 

At time t = 0, there is a singularity that can be incorporated in the above 
expressions. For simplicity, we only consider the expressions for t > 0. 

Since tsunameters have one component that measures the pressure at the 
bottom (see for example [H]), it is interesting to provide as well the expres- 
sion Pb{x, y, t) for the pressure at the bottom. The pressure p(x, y, z, t) can be 
obtained from Bernoulli's equation, which was written explicitly for the free 
surface in equation but is valid everywhere in the fluid: 

|^+l|V*P + (28) 

After linearization, equation ([28]) becomes 

Along the bottom, it reduces to 

^+g(-/i + C) + - = o, z = ~h. (30) 

ot p 

The time-derivative of the velocity potential is readily available in Fourier space. 
Inverting the Fourier and Laplace transforms and evaluating the resulting ex- 
pression &t z ~ —h gives for an instantaneous seabed deformation 



dt 



COS Lut ak. 



(27r)2 J J cosh^imh) 



The bottom pressure deviation from the hydrostatic pressure is then given by 

d(l) 



Pb{x,y,t) = - p-^^ 



- pgC- 

z=-h 



Away from the deformed seabed, ^ goes to so that pb simply is — p4>t\z=-h- 
The only difference between pb and pgrj is the presence of an additional cosh(TOft,) 
term in the denominator of pb- 
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3.2 Passive generation 

In this case equation pT|) becomes 

1^ = 0, z = -/i, (31) 
and the initial condition for t] now reads 

vix,y,0) = Ca{x,y), 

where (Q{x,y) is the seafloor deformation. Initial velocities are assumed to be 
zero. 

Again we apply the Fourier transform in the horizontal coordinates {x,y). 
The Laplace transform is not applied since there is no substantial dynamics in 
the problem. Equations ([9]), (|3ip and ([13]) become 

+ ^2)^=0, (32) 

^{kJ,-h,t)=Q, (33) 

^(fc,£,0,i)+ff|^(fc,^,0,t) = 0. (34) 
A general solution to Laplace's equation (|32|) is again given by 

(j){k, £, z, t) = A{k, I, t) cosh(mz) + B{k, I, t) sinh(mz), (35) 



where m = ^/ k"^ + P . The relationship between the functions A{k,£,t) and 
B{k, £, t) can be easily foimd from the boundary condition 



B{k, e, t) = A(k, e, t) tanh(m/i). (36) 
From equation (|34p and the initial conditions one finds 

A{k,e,t) ^ --Co{k,£)smu;t. (37) 
Substituting the expressions for the functions A and B in (|35p yields 

(j){k,£,z,t) = — — ^o(fc, ^) sinti>i^cosh(TOz) + tanh(m/i) sinh(mz)^ . (38) 
From p^. the free-surface elevation becomes 

r]{k,£,t) = (o{k,£) cosujt. 
Inverting the Fourier transform provides the general integral solution 

r]{x,y,t) = J j Coik,£)cosujt e'C^^+'^^^dfcdl (39) 
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Let us now evaluate the velocity field in the fluid domain. Equation 
gives the Fourier transform of the velocity potential (f>(x,y,z,t). Taking the 
Fourier transform of 



yields 



u{k,£,t;P) = -i(f>{kj,f3h,t)k 



= i—Coik, i) smu;tycosh{(3mh) + tanh(?Ti/i) sinh(/3m/i) j k. 

Inverting the Fourier transform gives the general formula for the horizontal 
velocities 

u{x,y,t;(3) = ^ ff kCo(fc,£)^^^(cosh(/?m/i)+tanh(?7i/i) sinh(/3TO/i))e*('=^+^y)dk 



Along the free surface f3 — 0, the horizontal velocity vector becomes 
u(x,y,t;0) = ^ //kCo(fc,^)^e'('^-+^-)dk. 



(40) 



Next we determine the vertical component of the velocity w{x,y,t; P) at a 
given vertical level z = (ih. It is easy to obtain the Fourier transform w{k, £, t; (3) 
by differentiating 



90 
dz 



sin Lot 



-mg- 



Co(fc, () ^sinh(/3TOft,)+tanh(TO/i) cosh(/3m/i) 



z=ph 

Inverting this transform yields 

m sin cut 



w{x,y,t;f3) = 



Caik,£)(smh{(3mh)+ 



tanhimh) cosh(/3m/i)) e^^'^^+^^'^dk 



for — 1 < /? < 0. Using the dispersion relation, one can write the vertical 
component of the velocity along the free surface {(3 = 0) as 



w{x,y,t;0) = 1 1 wBinwtCo{k,£)e^^'=+'^y^d\i. 



(41) 



All the formulas obtained in this section are valid only if the integrals converge. 
Again, one can compute the bottom pressure. At 2: = —h, one has 



dt 



9 ffCaik,e)e 



z— — h 



(2^)2 



i{kx+ly) 



cosh(m/i) 



■ cos Lot dk. 
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The bottom pressure deviation from the hydrostatic pressure is then given by 



Again, away from the deformed seabed, pb reduces to — p4>t\z=-h- The only 
difference between ph and pgrj is the presence of an additional cosh(m/i) term 
in the denominator of pb ■ 

The main differences between passive and active generation processes are 
that (i) the wave amplitudes and velocities obtained with the instantly moving 
bottom are lower than those generated by initial translation of the bottom 
motion and that (ii) the water column plays the role of a low-pass filter (compare 
equations ((25|) - (p7)) with equations p9|) - (|4T|) V High frequencies are attenuated 
in the moving bottom solution. Ward [36], who studied landslide tsunamis, also 
commented on the 1 / cosh(m/i) term, which low-pass filters the source spectrum. 
So the filter favors long waves. In the discussion section, we will come back to 
the differences between passive generation and active generation. 

3.3 Linear numerical method 

All the expressions derived from linear theory are explicit but they must be com- 
puted numerically. It is not a trivial task because of the oscillatory behaviour 
of the integrand functions. All integrals were computed with Filon type numer- 
ical integration formulas [6], which explicitly take into account this oscillatory 
behaviour. Numerical results will be shown in Section 6. 

4 Nonlinear shallow water equations 

Synolakis and Bernard [28j introduced a clear distinction between the various 
shallow-water models. At the lowest order of approximation, one obtains the 
linear shallow water wave equation. The next level of approximation provides 
the nondispersive nonlinear shallow water equations (NSW). In the next level, 
dispersive terms are added and the resulting equations constitute the Boussinesq 
equations. Since there are many different ways to go to this level of approxi- 
mation, there are a lot of different types of Boussinesq equations. The NSW 
equations arc the most commonly used equations for tsunami propagation (see 
in particular the code MOST developed by the National Oceanic and Atmo- 
spheric Administration in the US [29] or the code TUNAMI developed by the 
Disaster Control Research Center in Japan). They are also used for generation 
and runup/inundation. For wave runup, the effects of bottom friction become 
important and must be included in the codes. Our analysis will focus on the 
NSW equations. For simplicity, we assume below that h is constant. Therefore 
one can take h as reference depth, so that the seafloor is given by z = — 1 + e^. 



Pb{x,y,t) 



d(t> 



pgC- 



: = -h 
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4.1 Mathematical model 

In this subsection, partial derivatives are denoted by subscripts. When /x^ is a 
small parameter, the water is considered to be shallow. For the shallow water 
theory, one formally expands the potential 4> in powers of ji^: 

(/) = 00 + + + • • • . 

This expansion is substituted into the governing equation and the boundary 
conditions. The lowest-order term in Laplace's equation is 

00.. = 0. (42) 

The boundary conditions imply that 0o = 4'o{x^ t)- Thus the vertical velocity 
component is zero and the horizontal velocity components are independent of 
the vertical coordinate z at lowest order. Let us introduce the notation u := 
(f'Oxix, y,t) and v :~ (f>Qy(x,y,t). Solving Laplace's equation and taking into 
account the bottom kinematic condition yield the following expressions for 01 
and 02 : 

4ii{x,y,z,t) = -^Z'^{u.x+Vy) + z[Ct+e{uCx,+vCy)], 

1 / 1 \ 

02(x,y,z,t) = —z\Aux + Avy) + sie—\VCf~-Z''ACjiux + Vy) 

3 

-^Z^VC • V(m, + vy) - y (AC* + £A(wCx + vCy)) + 

z{-l + eC) [eVC • V(Ct + eiuCx + vCy)) - e'lVCpK + Vy)- 

^{-l + eO{ACt+£A{uC. + vCy)) , 

where 

Z=l + z-eC. 

The next step consists in retaining terms of requested order in the free- 
surface boundary conditions. Powers of e will appear when expanding in Taylor 
series the free-surface conditions around z = 0. For example, if one keeps terms 
of order e/i^ and /z^ in the dynamic boundary condition ([6]) and in the kinematic 
boundary condition ([4]), one obtains 

M^0ot - ^/(wte + Vty) + M^'? + ^£^^^{u^ + v^) = 0, (43) 

+ £{ur]x + vrjy) + (l + - C)) ("x + Vy) - (t - £{uC,x + vCy)] = 

\li\Aux + Avy). (44) 
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Differentiating first with respect to x and tlien with respect to y gives a 
set of two equations: 

Ut + e{uu., + VV^) +T]x- ^fJ.'^iutxx + Vtxy) = 0, (45) 

vt+ e{uuy + vvy) + T]y -]^^i^{utxy + vtyy) = 0. (46) 
The kinematic condition (|^^ becomes 

(77 - Ot + [^(1 + e(77 - C))]. + b(l + - C))]y = ^/i'(Au. + At;,). (47) 

Equations (|45p . (|^5)) and (|T7)) contain in fact various shallow- water models. The 
so-called fundamental NSW equations which contain no dispersive effects are 
obtained by neglecting the terms of order jj? : 

Ut + e{uux + vuy) +rix = 0, (48) 
vt + eiuvx + vvy)+ ijy = 0, (49) 
ilt + [u{l + e{v-C))]x + Hl + e{v-0)]y = Cf (50) 

Going back to a bathymetry h* {x, y, t) equal to 1 — eC,{x, y, t) and using the fact 
that (m, v) is the horizontal gradient of <j>o, one can rewrite the system of NSW 
equations as 

Ut + ^{u^ + v^)x + Vx = 0, (51) 
Vt + ^{u^ + v^)y + f]y = 0, (52) 

j^t + [uih* + srj)]x + Hh* + erj)]y = -h;. (53) 

The system of equations ([CT|) - ([55)) has been used for example by Titov and 
Synolakis for the numerical computation of tidal wave run-up. Note that 
this model does not include any bottom friction terms. 

The NSW equations with dispersion (|45|) -(|47 |) . also known as the Boussinesq 
equations, can be written in the following form: 

ut + ^iu"^ +v^)x+Vx'lf^^Aut = 0, (54) 

vt + ^{u^ +v\+ny-^fi^Avt = 0, (55) 

r]t + [u{h* +eT])]x + [v{h* +eT])]y-^fi^{Aux + Avy) = -^h;. (56) 

Kulikov ct al. [21] have argued that the satellite altimctry observations of the 
Indian Ocean tsunami show some dispersive effects. However the steepness 
is so small that the origin of these effects is questionable. Guesmia et al. [Tl] 
compared Boussinesq and shallow-water models and came to the conclusion that 
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the effects of frequency dispersion are minor. As pointed out in [TO], dispersive 
effects are necessary only when examining steep gravity waves, which are not 
encountered in the context of tsunami hydrodynamics in deep water. However 
they can be encountered in experiments such as those of Hammack |15j because 
the parameter ^ is much bigger. 



4.2 Numerical method 

In order to solve the NSW equations, a finite- volume approach is used. For 
example LeVeque |23] used a high-order finite volume scheme to solve a system 
of NSW equations. Here the flux scheme we use is the characteristic flux scheme, 
which was introduced by Ghidaglia et al. [9]. This numerical method satisfies 
the conservative properties at the discrete level. The NSW equations ([5T|) "([53 |) 
can be rewritten in the following conservative form: 

dw dF(w) dG(w) ^, 



dt dx dy 



where 



G 



{V,u,v), (58) 
+er,),|(u2 + t;2)+^,0) , (59) 

(^vih*+£^j),0,^{u^+v^)+Tj), (60) 

S = (-0,0,0). (61) 

The scheme we use is multi-dimensional by construction and does not re- 
quire the solution of any Riemann problem. For the sake of simplicity in the 
description of the numerical method, we assume that there is no y-variation 
and no source term S. Let us then consider a system of m-conservation laws 
(m > 1) 

dw dF(w) ^ , , 

- + ^=0, xeM, t>0, (62) 

where w £ R™ and F : M™ M™. Wc denote by A{w) the Jacobian matrix of 
F(w): 

dF 

Ay(w) = Tr^(w), l<i,j<m. (63) 



The system (|62p is assumed to be hyperbolic. In other words, for every w 
there exists a smooth basis (ri(w), . . . , rm(w)) of consisting of eigenvectors 
of A(w). Said differently, there exists Afe(w) e M such that A(w)rfc(w) = 
Afc(w)rfc(w). It is then possible to construct (£i(w), . . . ,£„i{w)) such that 

*A(w)4(w) = Afc(vir)4(w) and 4(w) • rp(w) Skp- 

Let M = Uj£z[xj-i/2, Xj+1/2] be a one-dimensional mesh. Let also M+ = 
U„gM[in, ^Ti+i]- Let us discretize (j62p by a finite volume method. We set 



4.2 Numerical method 



21 



and 



The system ((62|) can then be rewritten (exactly) as 

At 



(F;Vi/2 - F^i/2)- (64) 



For a three-point cxphcit numerical scheme one has 

F}Vi/2«f;Xw;,w^i), (65) 
where the function f is to be specified. Multiplying ([62)) by A(w) yields 

This shows that the flux F(w) is advected by A(w). The numerical flux 
f"(w", w"_|_]^) represents the flux at an interface. Using a mean value /^"_|_i/2 of 
w at this interface, we replace by the linearization 

^IM + A(n- .)^-0 (67) 

We define the fc— th characteristic flux component to be Fk{w) — ^fc(/^"+i/2) ' 
F(w). It follows that 

+ ^'=(/^.+i/2)^^ - 0. (68) 

This linear equation can be solved explicitly for Ff.(w). As a result it is nat- 
ural to dcflne the characteristic flux f-"^ at the interface between two cells 
[xj_i/2, a;j+i/2] and [xj+1/2, Xj+3/2] as follows: for fc e 1, . . . , to, 



■ ff^'"(w;',w;'+i) = 4(mJ+i/2) • F(w7), when Xki^i]+,/2) > 0, 



WU1/2) ■ ff^'"(w;^w;Vl) = WUir2) ■ f(w;Vi), when Afc(^;Vi/2) < o, 

«fe(Aij + l/2) • Ij ' Wj + l) - «fe(Aij + l/2) • 1 ^ I ) 



when Afe(/x^*_^_^/2) = 0- Here 



The characteristic flux can be written as 



f--'"(w;,w;Vi) = f^"(A^";w7,w;Vi) 
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where 

fCFt , F(wi)+F(w2) F(w2)-F(wi) 

f (/^;wi,W2) = sgn(A(/i(wi,W2)) . (69) 

The sign of the matrix A(/i) is defined by 

k—ni 

sgn(A(M))$ = ^ sgn(Afc)(4(/^) • $)rfc(M). 



Going back to (|64p . one can construct the following explicit scheme: 

= - ^ (ff^-(w;, w;V,) - ff ^-(w;^,, w^)) . (70) 

The characteristic flux scheme (j69p gives very good results, especially when 
complex systems are considered [3]. In our case, we have to consider equation 
(|62[) in two dimensions and to discretise the source term too: 

One can refer to [TU] for these two extensions. 



5 Numerical method for the full equations 

The fully nonlinear potential flow (FNPF) equations Q-® are solved by using 
a numerical model based on the Boundary Element Method (BEM) . An accurate 
code was developed by Grilli et al. [H]. It uses a high-order three-dimensional 
boundary element method combined with mixed Eulerian-Lagrangian time up- 
dating, based on second-order explicit Taylor expansions with adaptive time 
steps. The efficiency of the code was recently greatly improved by introducing 
a more efficient spatial solver, based on the fast multipole algorithm [7]. By 
replacing every matrix-vector product of the iterative solver and avoiding the 
building of the influence matrix, this algorithm reduces the computing complex- 
ity from 0{N^) to nearly 0{N) up to logarithms, where N is the number of 
nodes on the boundary. 

By using Green's second identity, Laplace's equation ([T|) is transformed into 
the boundary integral equation 

a(xO0(x,) = ^ (^|^(x)G(x,xO - 0(x)|^(x,xo) dT, (72) 

where G is the three-dimensional free space Green's function. The notation 
dG/dn means the normal derivative, that is dG/dn = VG • n, with n the 
unit outward normal vector. The vectors x = {x,y,z) and x; ~ {xi,yi,zi) are 
position vectors for points on the boundary, and Q!(x() = 0(/(47r) is a geometric 
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coefficient, with 0i the exterior sohd angle made by the boundary at point X;. 
The boundary T is divided into various parts with difi'erent boundary conditions. 
On the free surface, one rewrites the nonhncar kinematic and dynamic boundary 
conditions in a mixed Eulerian-Lagrangian form, 

D(h \ 

= -.9Z+-V0-V0, (74) 

with R the position vector of a free-surface fluid particle. The material deriva- 
tive is defined as 

For time integration, second-order explicit Taylor series expansions are used 
to find the new position and the potential on the free surface at time t + At. 
This time stepping scheme presents the advantage of being explicit, and the 
use of spatial derivatives along the free surface provides a better stability of the 
computed solution. 

The integral equations are solved by BEM. The boundary is discrctizcd into 
TV collocation nodes and M high-order elements are used to interpolate between 
these nodes. Within each clement, the boundary geometry and the field vari- 
ables arc discrctizcd using polynomial shape functions. The integrals on the 
boundary are converted into a sum on the elements, each one being calculated 
on the reference element. The matrices are built with the numerical computa- 
tion of the integrals on the reference element. The linear systems resulting from 
the two boundary integral equations (one for the pair (0, d(f)/dn) and one for 
the pair [dcj)/ dt,d^4>/ dtdn)) are full and non symmetric. Assembling the ma- 
trix as well as performing the integrations accurately arc time consuming tasks. 
They are done only once at each time step, since the same matrix is used for 
both systems. Solving the linear system is another time consuming task. Even 
with the GMRES algorithm with preconditioning, the computational complex- 
ity is 0{N^)., which is the same as the complexity of the assembling phase. The 
introduction of the fast multipolc algorithm reduces considerably the complex- 
ity of the problem. The matrix is no longer built. Far away nodes are placed 
in groups, so less time is spent in numerical integrations and memory require- 
ments are reduced. The hierarchical structure involved in the algorithm gives 
automatically the distance criteria for adaptive integrations. 

Grilli et al. [13] used the earlier version of the code to study tsunami gen- 
eration by underwater landslides. They included the bottom motion due to the 
landslide. For the comparisons shown below, wc only used the passive approach: 
we did not include the dynamics of the bottom motion. 



6 Comparisons and discussion 

The passive generation approach is followed for the numerical comparisons be- 
tween the three models: (i) linear equations, (ii) NSW equations and (iii) fully 
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nonlinear equations. As shown in Section 3, this generation process gives the 
largest transient-wave amplitudes for a given permanent deformation of the 
seafloor. Therefore it is in some sense a worst case scenario. 

The small dimensionless numbers £ and //^ introduced in ([7]) represent the 
magnitude of the nonlinear terms and dispersive terms in the governing equa- 
tions, respectively. Hence, the relative importance of the nonlinear and the 
dispersive effects is given by the parameter 

nonlinear terms e a}? ,_ ^ 

^ = T = — = O's 

dispersive terms /i^ dfi 

which is called the Stokes (or Ursell) number [34] An important assumption 
in the derivation of the Boussinesq system (j54p - ((56|) is that the Ursell number is 
0(1). Here, the symbol O(-) is used informally in the way that is common in the 
construction and formal analysis of model equations for physical phenomena. 
We are concerned with the limits e — > and /x ^ 0. Thus, S = 0(1) means that, 
as e ^ and jj. ^ 0, S takes values that are neither very large nor very small. 
We emphasize here that the Ursell number does not convey any information 
by itself about the separate negligibility of nonlinear and frequency dispersion 
effects. Another important aspect of models is the time scale of their validity. 
In the NSW equations, terms of order 0{e^) and 0{fi^) have been neglected. 
Therefore one expects these terms to make an order-one relative contribution 
on a time scale of order min(£~^, 

All the figures shown below are two-dimensional plots for convenience but we 
recall that all computations for the three models are three-dimensional. Figure 
[6] shows profiles of the free-surface elevation along the main direction of propaga- 
tion (y— axis) of transient waves generated by a permanent seafloor deformation 
corresponding to the parameters given in Table [1] This deformation, which has 
been plotted in figure [21 has been translated to the free surface. The water 
depth is 100 m. The small dimensionless numbers are roughly e = 5xl0~^ and 
/X = 10~^, with a corresponding Ursell number equal to 5. One can see that the 
front system splits in two and propagates in both directions, with a leading wave 
of depression to the left and a leading wave of elevation to the right, in qual- 
itative agreement with the satellite and tide gauge measurements of the 2004 
Sumatra event. When tsunamis are generated along subduction zones, they 
usually split in two; one moves quickly inland while the second heads toward 
the open ocean. The three models arc almost undistinguishable at all times: the 
waves propagate with the same speed and the same profile. Nonlinear effects 
and dispersive effects are clearly negligible during the first moments of transient 

^One finds sometimes in the literature a subtle difference between the Stokes and Ursell 
numbers. Both involve a wave amplitude multiplied by the square of a wavelength divided by 
the cube of a water depth. The Stokes number is defined specifically for the excitation of a 
closed basin while the Ursell number is used in a more general context to describe the evolution 
of a long wave system. Therefore only the characteristic length is different. For the Stokes 
number the length is the usual wavelength A related to the frequency id hy X 2TTy/~gd/u}. 
In the Ursell number, the length refers to the local wave shape independent of the exciting 
conditions. 
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Wave profile at time t = s. Slice of 2D surface. 




1,1 \ I ^ I I I I I 
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Figure 6: Comparisons of the free-surface elevation at x = resulting from the 

integration of the linear equations (•••), NSW equations ( ) and nonlinear 

equations (— ) at different times of the propagation of transient waves generated 
by an earthquake {t = s, t = 95 s, t ~ 143 s, t ~ 191 s). The parameters for 
the earthquake are those given in Table 1. The water depth is h = 100 m. One 
has the following estimates: e = 5 x 10^"*, fi^ = 10^^ and consequently 5 = 5. 
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Wave profile at time t = 52 s 




y (the msin diret^lion of propagation, in l<m) 



Figure 7: Comparisons of the free-surface elevation at x = resulting from the 

integration of the linear equations (■■■), NSW equations ( ) and nonlinear 

equations (— ) at different times of the propagation of transient waves generated 
by an earthquake {t = 52 s, t = 104 s, t = 157 s). The parameters for the 
earthquake are those given in Table 1. The water depth is h = 500 m. One has 
the following estimates: e = 10^^, /i^ = 2.5 x 10^'^ and consequently S = 0.04. 

waves generated by a moving bottom, at least for these particular choices of e 
and n- 

Let us now decrease the Ursell number by increasing the water depth. Figure 
[7| illustrates the evolution of transient water waves computed with the three 
models for the same parameters as those of figure El except for the water depth 
now equal to 500 m. The small dimensionless numbers are roughly e = 10^'^ 
and fi = 5x 10~^, with a corresponding Ursell number equal to 0.04. The linear 
and nonlinear profiles cannot be distinguished within graphical accuracy. Only 
the NSW profile is slightly different. 

Let us introduce several sensors (tide gauges) at selected locations which are 
representative of the initial deformation of the free surface (see figure |8]). One 
can study the evolution of the surface elevation during the generation time at 
each gauge. Figure ^ shows free-surface elevations corresponding to the linear 
and nonlinear shallow water models. They are plotted on the same graph for 
comparison purposes. Again there is a slight difference between the linear and 
the NSW models, but dispersion effects are still small. 

Let us decrease the Ursell number even further by increasing the water depth. 
Figures [10] and [Tl] illustrate the evolution of transient water waves computed 
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Figure 8: Top view of the initial free surface deformation showing the location 
of six selected gauges, with the following coordinates (in km): (1) 0,0 ; (2) 0,3 ; 
(3) 0,-3 ; (4) 10,5; (5) —2,5 ; (6) 1,10. The lower oval area represents the initial 
subsidence while the upper oval area represents the initial uplift. 



6 Comparisons and discussion 28 

























Tide gaug 





































































































20 40 



100 120 




20 




Figure 9: Transient waves generated by an underwater earthquake. Compar- 
isons of the free-surface elevation as a function of time at the selected gauges 

shown in figure [S] — , linear model ; nonlinear shallow water model. The 

time t is expressed in seconds. The physical parameters are those of figure [7l 
Since the fully nonlinear results cannot be distinguished from the linear ones, 
they are not shown. 
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with the three models for the same parameters as those of Figure El except 
for the water depth now equal to 1 km. The small dimensionless numbers 
are roughly e = 5 x 10^^ and /.j = 0.1, with a corresponding Ursell number 
equal to 0.005. On one hand, linear and fully nonlinear models are essentially 
undistinguishable at all times: the waves propagate with the same speed and the 
same profile. Nonlinear effects are clearly negligible during the first moments 
of transient waves generated by a moving bottom, at least in this context. On 
the other hand, the numerical solution obtained with the NSW model gives 
slightly different results. Waves computed with this model do not propagate 
with the same speed and have different amplitudes compared to those obtained 
with the linear and fully nonlinear models. Dispersive effects come into the 
picture essentially because the waves are shorter compared to the water depth. 
As shown in the previous examples, dispersive effects do not play a role for long 
enough waves. 

Figure [12] shows the transient waves at the gauges selected in figure [8] One 
can see that the elevations obtained with the linear and fully nonlinear models 
are very close within graphical accuracy. On the contrary, the nonlinear shallow 
water model leads to a higher speed and the difference is obvious for the points 
away from the generation zone. 

These results show that one cannot neglect the dispersive effects any longer. 
The NSW equations, which contain no dispersive effects, lead to different speed 
and amplitudes. Moreover, the oscillatory behaviour just behind the two front 
waves is no longer present. This oscillatory behaviour has been observed for 
the water waves computed with the linear and fully nonlinear models and is 
due to the presence of frequency dispersion. So, one should replace the NSW 
equations with Boussinesq models which combine the two fundamentals effects 
of nonlincarity and dispersion. Wei et al. |37| provided comparisons for two- 
dimensional waves resulting from the integration of a Boussinesq model and 
the two-dimensional version of the FNFF model described above. In fact they 
used a fully nonlinear variant of the Boussinesq model, which predicts wave 
heights, phase speeds and particle kinematics more accurately than the standard 
weakly nonlinear approximation first derived by Peregrine [27j and improved 
by Nwogu's modified Boussinesq model [25]. We refer to the review [20] on 
Boussinesq models and their applications for a complete description of modern 
Boussinesq theory. 

From a physical point of view, we emphasize that the wavelength of the 
tsunami waves is directly related to the mechanism of generation and to the 
dimensions of the source event. And so is the dimensionless number fi which 
determines the importance of the dispersive effects. In general it will remain 
small. 

Adapting the discussion by Bona ct al. [2], one can expect the solutions 
to the long wave models to be good approximations of the solutions to the 
full water-wave equations on a time scale of the order min(£~^, and also 
the neglected effects to make an order-one relative contribution on a time scale 
of order min(e~^, e^^ju"^). Even though we have not computed precisely 
the constant in front of these estimates, the results shown in this paper are in 
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Figure 10: Comparisons of the free-surface elevation at x = resulting from 

the integration of the linear equations (— • — ), NSW equations ( ) and FNPF 

equations (— ) at different times of the propagation of transient waves generated 
by an earthquake (t = 50 s, i = 100 s). The parameters for the earthquake 
are those given in Table 1. The water depth is 1 km. One has the following 
estimates: e = 5 x 10^^, fj,'^ = 10~^ and consequently S = 0.005. 
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Figure 11: Same as figure fTOl for later times (t = 150 s, i = 200 s). 
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Figure 12: Transient waves generated by an underwater earthquake. The phys- 
ical parameters are those of figures [10] and [11] Comparisons of the free-surface 
elevation as a function of time at the selected gauges shown in figure]!] — , linear 

model ; nonlinear shallow water model. The time t is expressed in seconds. 

The FNPF results cannot be distinguished from the linear results. 
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agreement with these estimates. Considering the 2004 Boxing Day tsunami, 
it is clear that dispersive and nonhnear effects did not have sufficient time to 
develop during the first hours due to the extreme smallness of e and /i^, except 
of course when the tsunami waves approached the coast. 

Let us conclude this section with a discussion on the generation methods, 
which extends the results given in [5] H. We show the major differences between 
the classical passive approach and the active approach of wave generation by a 
moving bottom. Recall that the classical approach consists in translating the 
sea bed deformation to the free surface and letting it propagate. Results are 
presented for waves computed with the linear model. 

Figure [T2] shows the waves measured at several artificial gauges. The pa- 
rameters are those of Table [U and the water depth is = 500 m. The solid 
line represents the solution with an instantaneous bottom deformation while the 
dashed line represents the passive wave generation scenario. Both scenarios give 
roughly the same wave profiles. Let us now consider a slightly different set of 
parameters: the only difference is the water depth which is now h = 1 km. As 
shown in figure fT4], the two generation models differ. The passive mechanism 
gives higher wave amplitudes. 

Let us quantify this difference by considering the relative difference between 
the two mechanisms defined by 



Intuitively this quantity represents the deviation of the passive solution from 
the active one with a moving bottom in units of the maximum amplitude of 
Va.ctive{x, y,t). 

Results are presented on figures (fT5|) and (fT6|) . The differences can be easily 
explained by looking at the analytical formulas (|25p and ([39|) of Section 3. These 
differences, which can be crucial for accurate tsunami modelling, are twofold. 

First of all, the wave amplitudes obtained with the instantly moving bottom 
are lower than those generated by the passive approach (this statement follows 
from the inequality coshmft, > I). The numerical experiments show that this 
difference is about 6% in the first case and 20% in the second case. 

The second feature is more subtle. The water column has an effect of a low- 
pass filter. In other words, if the initial deformation contains high frequencies, 
they will be attenuated in the moving bottom solution because of the presence of 
the hyperbolic cosine cosh(m/i) in the denominator which grows exponentially 
with TO. Incidcntly, in the framework of the NSW equations, there is no differ- 
ence between the passive and the active approach for an instantaneous seabed 
deformation [5^155] . 

If we prescribe a more realistic bottom motion as in [3] for example, the 
results will depend on the characteristic time of the seabed deformation. When 
the characteristic time of the bottom motion decreases, the linearized solution 

^In figures 1 and 2 of [5], a mistake was introduced in the time scale. All times must be 
multiplied by a factor VlOOO. 
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Figure 13: Transient waves generated by an underwater earthquake. The com- 
putations are based on Unear wave theory. Comparisons of the free-surface 
elevation as a function of time at selected gauges for active and passive genera- 
tion processes. The time t is expressed in seconds. The physical parameters are 
those of figure [71 In particular, the water depth is ft- = 500 m. 
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Figure 14: Same as figure [T31 except for tiie water deptii, wiiicii is equal to 1 
km. 

tends to the instantaneous wave generation scenario. So, in the framework of 
linear water wave equations, one cannot exceed the passive generation amplitude 
with an active process. However, during slow events, Todorovska and Trifunac 
|31| have shown that amplification of one order of magnitude may occur when the 
sea fioor uplift spreads with velocity similar to the long wave tsunami velocity. 

7 Conclusions 

Comparisons between linear and nonlinear models for tsunami generation by an 
underwater earthquake have been presented. There are two main conclusions 
that are of great importance for modelling the first instants of a tsunami and for 
providing an efficient initial condition to propagation models. To begin with, a 
very good agreement is observed from the superposition of plots of wave profiles 
computed with the linear and fully nonlinear models. Secondly, the nonlinear 
shallow water model was not sufficient to model some of the waves generated 
by a moving bottom because of the presence of frequency dispersion. However 
classical tsunami waves are much longer, compared to the water depth, than the 
waves considered in the present paper, so that the NSW model is also sufficient 
to describe tsunami generation by a moving bottom. Comparisons between the 
NSW equations and the FNPF equations for modeling tsunami run-up are left 
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Figure 15: Relative difference between the two solutions shown in figure 1131 
The time t is expressed in seconds. 




Figure 16: Relative difference between the two solutions shown in figure [Til 
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for future work. Another aspect which deserves attention is the consideration of 
Earth rotation and the derivation of Boussinesq models in spherical coordinates. 
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